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Abstract. We study the regular or chaotic nature of motion in a disk galaxy 
with a dense nucleus and an asymmetric dark halo. Two cases, the 2D model 
and the 3D model, are investigated. In the 2D model, a considerable fraction of 
the phase plane is covered by chaotic orbits. Two factors seem to be responsible 
for the chaotic motion: (i) the dense nucleus and (ii) the asymmetries in the 
dark halo. Our numerical experiments suggest, that there are several chaotic 
components on the Poincarc phase plane. Different chaotic components arc 
induced by the asymmetries in the halo. Each chaotic component seems to 
have a different value of the Lyapunov Characteristic Exponent, for small values 
of the asymmetry parameter A and a unique LCE for larger values of A. A 
comparison of the present results with outcomes from previous work is also 
presented. 
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1. INTRODUCTION 

Galaxies are often surrounded by halos. In most cases the shape of halo is 
spherical or nearly spherical but there are also indications that the shape of halo 
may be a biaxial or even a triaxial ellipsoid (see Kunihito et al. 2000; Oiling & 
Merrifield 2000; Wechsler et al. 2002; Caranicolas & Zotos 2009). The distinction 
between the halo and the main body of the galaxy is clearest in disk galaxies, 
where the spherical or the triaxial shape of the halo contrasts with the flat disc. 
In an elliptical galaxy, there is no sharp transition between the body of the galaxy 
and the halo. The visible part of the halo, is occupied by population II objects, 
including globular clusters and old individual stars. Beyond this, is a much larger 
region, called the dark halo or extended halo, containing large amounts of dark 
matter. Note that the stellar halo is not an inner part of the dark halo. They 
are two physically distinct components with different formation histories. We 
would like to point out that the above mentioned triaxial dark halo models are all 
symmetrical, unlike the dark halo model in the present paper. 

The need for existence of dark halos results from many independent sources 
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such as rotation curves, hot gas in clusters, velocity dispersions of galactic groups, 

gravitational lensing and microwave radiation. The presence of dark matter in 
the halo is demonstrated by its gravitational effect on the rotation curve of the 
disk. Without large amounts of mass in the extended halo, the rotational velocity 
of the galaxy should decrease at large distance from the galactic core. However, 
observations of disk galaxies, particularly radio observations of line emission from 
neutral atomic hydrogen, show that the rotation curve of most spiral galaxies 
remains flat far beyond the visible matter. The absence of any visible matter to 
account for these observations implies the presence of unobserved matter. The 
nature of dark matter in the galactic halo of spiral galaxies is still undetermined, 
but the most popular theories accept that the dark halo is populated by vast 
numbers of small bodies known as MACHOs. Observations of the halo of the 
Milky Way show that the number of MACHOs is not likely to be sufficient to 
account for the required mass. 

No doubt that the behavior of orbits in the galactic disk is affected by the 
presence of the halo. Of particular interest is to study the effect of the asymmetries 
of the dark halo in the behavior of orbits in the disk. In order to investigate this 
behavior we have constructed a composite 3D dynamical model of a disk galaxy. 
The model is described in Section 2. In Section 3 we use the x — Px phase plane 
in order to study the motion in the 2D model. In Section 4 we use the results 
obtained from the 2D model in order to visualize the motion in the 3D model. Of 
special interest is to find the behavior of 3D orbits in different chaotic components 
observed in the 2D model. We close our research with discussion and conclusions 
which are presented in Section 5. 



2. PRESENTATION OF THE DYNAMICAL MODEL 

Our model consists of three distinct components - the disk halo, the nucleus 
and the dark halo component. The disk halo component is described by the 
potential 

Vd{x,y,z)= (1) 

J }? + x"^ + y"^ + (a + V/i^ + ^2)2 

where Md is the mass, h is the core radius of the disk halo, a is the disk scale 
length and h is the disk scale height. Potential (1) was introduced by Miyamoto 
& Nagai (1975). The dense massive nucleus is described by the spherical potential 

Vn{x,y,z) = —==^^===, (2) 
^^ x^ + y-^ + z-^ + cf^ 

where M„ is the mass and c„ is the scale length of the nucleus. For the dark halo 
we use the logarithmic potential 

Vh{x, y, z) = "^Inlx^ +y^ + z^- \x^ + cl]. (3) 

Here Ch is the scale length of the dark halo component and the parameter vq is 

used for the consistency of the galactic units. The term —\x^,\ ^ 1 represents 
the deviation of the halo from spherical symmetry (see also Binney & Tremaine 
2008). 
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We use the system of galactic units, where the unit of length is 1 kpc, the unit of 
mass is 2.325 x lO"^ Mq and the unit of time is 0.97748 x 10^ yr. The velocity imit is 
10 km/s, while G is equal to unity. The energy unit (per unit mass) is 100 (km/s)^. 
In the above units we use the values: vq = 20, a = 3,b = 6, h = 0.325, Md = 12000. 
The value of c„ and Ch is 0.25 and 8 respectively while M„ and A are treated as 
parameters. 

The total potential responsible for the motion of a test particle of unit mass in 
the galaxy is 

VTix,y,z) = Vd + Vn + Vh. (4) 

There are three lines of arguments to justify the choice of potential (4). The 
first is that the disk nucleus potential system given by equations (1) and (2) is a 
classical realistic potential, which describes very well the motion in an active disk 
galaxy. The second is that potential (3) describes satisfactory the deviation of the 
dark halo from spherical symmetry, and the third is that there is observational 
evidence that asymmetries in the halos do exist (see Xu et al. 2007). However we 
must emphasize that (Xu et al. 2007) investigated the stellar halos, not the dark 
halos. 

The equations of motion are: 

..__dVT ..__dVT ..__dVT 
^ dx dy dz ' 

where the dot indicates derivative with respect to the time. The corresponding 
Hamiltonian is 

H=^ipl+pl+pl) + VTix,y,z)=E, (6) 

where Px^Py^Pz are the momenta per unit mass, conjugate to x, y and z while E 
is the numerical value of the Hamiltonian. 



3. NUMERICAL RESULTS FOR THE 2D MODEL 

In this case we integrate numerically the equations of motion (5) with z = Pz = 
0. The corresponding Hamiltonian is 

H2 = l{pl+pl) + VT{x,y)=h2, (7) 

where /i2 is the numerical value of the energy of the test particle. 

We believe that it is a good idea to start from the 2D system for two basic 
reasons: (i) in the 2D system we can use the x ^ px,y = 0,Py > Poincare phase 
plane in order to locate the areas of regular and chaotic motion, and (ii) we can 
use the experience gained from the study of the 2D system in order to explore a 
more complicated 3D system. 

Figure 1 shows the x — px phase plane when A = 0.015, M„ = 400, /i2 = 500. 
We see that the larger part of the phase plane is covered by regular orbits. There 
are also parts of the phase plane occupied by chaotic orbits. Three distinct areas 
of chaotic motion are observed: a considerable chaotic layer near the central region 
and two secondary chaotic components which appear near the unstable periodic 
orbit produced by secondary resonances. Figure 2 is similar to Figure 1 but for 
/i2 = 300. Here one can see only the chaotic layer near the central part of the 
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phase plane while no other chaotic component seems to be present. The rest of the 

phase plane is covered by regular orbits. Thus we may conclude that, as a result, 
the asymmetry in the dark halo has to produce a different chaotic component 
for the high energy stars. In other words, the asymmetry affects not only stars 
approaching the dense nucleus but also stars moving far from the nuclear region. 
On the other hand, from the totality of low energy stars only those approaching 
the nucleus arc on chaotic orbits. 

Figure 3 is similar to Figure 1 but for A = 0.03. In Figure 3 we see that 
the three chaotic components have merge to produce a large chaotic sea, while 
the small regular regions are confined near the stable resonant periodic points. 
Furthermore, in Figure 4 we see that the picture is almost the same as that shown 
in Figure 2. The only difference is that here one can observe some additional small 
islands produced by secondary resonances. Thus, one can see that the increase of 
the asymmetry in the halo affects drastically the high energy stars, producing large 
chaotic regions by merging different chaotic components. Moreover, low energy 
stars seem to display only one chaotic layer, which is produced by low energy stars 
approaching the dense nucleus. 

Figures 5(a-d) show four typical orbits for the 2D potential. The values of 
the parameters are as in Figure 1. Panel (a) shows a quasi periodic orbit. The 
initial conditions are: xg = — 6.5,j/o = Q,Px{) = 19, /?.2 = 500. The value of py 
is always found from energy integral (7). Panel (b) shows a tube orbit with the 
initial conditions xq = 7.0, j/o = 0, Pxo = and /12 = 300. A quasi periodic orbit 
is also presented in panel (c). The initial conditions are: = -3.0, yo = PxO = 
and /i2 = 300. Note that all the above quasi periodic orbits do not approach the 
dense nucleus. Panel (d) shows a chaotic orbit with the initial conditions xq = 
0.5, yo = PxO = and /12 = 500. This orbit comes arbitrary close to the nucleus. 
The orbits were calculated for a time period of 100 time units. 

It would be of particular interest to follow the evolution of chaotic regions, as 
the parameter A increases or as the mass of the nucleus M„ increases, when all 
other parameters are kept constant. The results are presented in Figures 6 and 7. 
Figure 6 shows a plot of the percentage of the surface of section A% covered by 
chaotic orbits vs. A, when M„ = 400. The values of all other parameters are as in 
Figure 1. Dots indicate the values found numerically, while the solid line is a two 
degree polynomial fit. We see that A% increases rapidly as A increases. Figure 7 
shows a plot of A% vs. M„, when A = 0.015. The values of all other parameters 
are as in Figure 1. Again, dots indicate the values found numerically, while the 
solid line is a two degree polynomial fit. Here we see that A% increases as M„ 
increases but at a smaller rate. 

It is also of interest to compute the Lyapunov Characteristic Exponents (LCE) 
(see Lichtenberg & Lieberman 1992) in order to estimate the degree of chaos in 
diflFerent chaotic regions of Figure 1, i.e., when A = 0.015, M„ = 400, /12 = 500. 
The LCE in all cases was computed for a time period of 20000 time imits. Note 
that this time period is at least 100 times larger than the corresponding Hubble 
time. The LCE was found to be 0.195 in the chaotic region near the center. 
The value of the LCE on the left and right chaotic components was found to 
be 0.020 and 0.010, respectively. As expected, each chaotic component has a 
diflferent LCE. This phenomenon was observed in disk galaxies withoiit a halo 
component (see Fig. 3 in Papadopoulos & Caranicolas 2005). It was also observed 
in potentials made up of perturbed harmonic oscillators (see Saito and Ichimura 
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Fig. 1. The x — Px phase plane of the 2D system. The values of the parameters are: 
uo = 20,a = 3,& = 6,^ = 0.325, Md = 12000, c„ = 0.25, Ch = 8,M„ = 400, A = 0.015 
and h2 = 500. 
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Fig. 2. The same as in Figure 1 but for h2 = 300. 
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Fig. 3. The x — px phase plane of the 2D system. The values of the parameters are: 
Wo = 20, a = 3, fe = 6, /i = 0.325, Md = 12000, c„ = 0.25, ch = 8, M„ = 400, A = 0.03 and 
/i2 = 500. 
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Fig. 4. The same as in Figure 3 but for h2 = 300. 
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Fig. 5. Panels (a)-(d): typical orbits for the 2D potential. The values of the param- 
eters are as in Figure 1. Panel (a) - upper left: a quasi-periodic orbit with the initial 
conditions: xq = —6.5, yo = 0, p^o = 19 and /i2 = 500. Panel (b) - upper right: a 
tube orbit with the initial conditions: xq = 7.0, yo = 0, Pxo = and h2 — 300. Panel 
(c) - lower left: a quasi periodic orbit with the initial conditions xq — —3.0, yo = 0, 
PxO = and /i2 ~ 300. Panel (d) - lower right: a chaotic orbit with the initial conditions 
xo = 0.5, yo ~ 0, PxO ~ and h2 = 500. The value of py is always found from the energy 
integral. 

1978). Moreover, the LCE for all orbits, starting in the chaotic sea of Figure 3, 
i.e., when A = 0.03, Af„ — 400, /12 — 500, was found to be 0.165. This result is not 
surprising as in Figure 3 different chaotic regions have merged. 

4. ORBITS IN THE 3D MODEL 

Let us now turn to the properties of orbits in the 3D potential. In order to do 
this, we compute orbits with initial conditions {xo,PxO, zq), yo = Pzo = 0, where 
{xo,Pxo) is a point inside the limiting curve of the 2D system. The limiting curve 
corresponds to 

^pI + Vt{x) = h2 (8) 

and can be obtained from (7) if we set y = Py = 0. For convenience, we take 
E — h2 and the value of Pyo is always found from the energy integral (6). 

Our numerical calculations suggest that all orbits in the 3D model, starting 
with initial conditions (xo^Pxo) on the chaotic zones of the 2D system of Figure 1, 
i.e., when A = 0.015, M„ = 400, = /12 = 500, and all permissible values of zq, 
i.e. such as to obtain real values for py, are chaotic. The values of the LCE are 
different for the 3D orbits starting in different chaotic zones. Thus, the value of 
the LCE for orbits starting with {xo,Pxo) in the chaotic zone near the center was 
found to be 0.105. For the chaotic zone on the left-hand side of Figure 1, the LCE 
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Fig. 6. A plot of the percentage of the surface of section AVo covered by cliaotic orbits 
vs. A, when Af„ = 400. 
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Fig. 7. A plot of the percentage of the surface of section A% covered by chaotic orbits 
vs. M„, when A = 0.015. 



was found to be 0.015 while for the chaotic zone on the right-hand side of Figure 
1 it was 0.008. These results seem to be correct because there is no evidence 
yet that in 3D systems with divided phase space a completely connected chaotic 
component actually exists (see Cincotta et al. 2006). On the other hand, all orbits 
with initial conditions {xo,PxQ, zq), j/q — Pzo — 0, with values of {xo,Pxo) in the 
chaotic sea of Figure 3, were found chaotic for all permissible values of zq. The 
corresponding LCE was found to converge to a common value, equal to 0.124. 

In order to have a complete picture of the behavior of orbits in the 3D system, 
we have calculated a large number of orbits with initial conditions, (xq^Pxo, zo)^ 
Vo = PzO = with the values of {xq^p^o) in the regular regions of Figure 1. It 
was found that for the values \zo\ < 4 all tested orbits were found regular, while 
when Izol > 4 the orbits were found chaotic. The value of the LCE was found to be 



Chaos in a disk galaxy model induced by dark halo 



213 




Fig. 8. Panels (a)-(d): representative orbits in the 3D potential. Initial conditions 
for xo,Pxo as in Figures 5a-5d, respectively, while zg — 0.1 and E — h2 for all orbits. 
Note that the chaotic orbit in Figure 5d is scattered to the halo. 



between 0.005-0.008 for orbits starting on the left-hand side of Figure 1, while 
for orbits on the right-hand side of Figure 1 the LCE takes values in a range of 
0.010-0.015. 

Figures 8 (a-d) show four representative orbits in the 3D potential. In order 
to better visualize the behavior of the 3D orbits, we have used the same initial 
conditions of the orbits as in Figures 5 (a-d) for the 2D system, with the same 
value of energy, E — h2- The value of zq was taken to be equal to 0.1 for all orbits. 
It is seen that all orbits retain their regular or chaotic character. Furthermore, 
one can see only chaotic orbits approach the dense nucleus while regular orbits do 
not. All orbits were calculated for a time period of 100 time units. 

It is also interesting to note that, from all the above four orbits starting near the 
galactic plane, only that approaching the nucleus is deflected to the halo. There- 
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Fig. 9. Evolution of the angular momentum vs. time: the left panel is for the quasi- 
periodic orbit of Figure 8b and the right panel is for the chaotic orbit of Figure 8d. 

fore, we can say that orbits moving near the galactic plane, when approaching the 
dense nucleus, are scattered to higher z. This behavior is similar to that observed 
about twenty years ago by Caranicolas & Inannen (1991) in axially symmetrical 
disk galaxies, hosting dense nuclei. Today it is well known that for this scattering is 
responsible the strong Fz nuclear force, that acts upon the low angular momentum 
stars (see also Caranicolas & Papadopoulos 2003). 

We would also like to remind the reader that, in any case, the orbits that 
are scattered to the halo are orbits of low angular momentum. As here the 
component of the angular momentum is not conserved, we can compute its mean 
value 

1 " 

<L, >=-VL,fc. (9) 
n ^ — ' 

Figure 9a shows the evolution of angular momentum vs. time for the orbit 
shown in Figure 8b. As expected, one can observe a quasi periodic curve while 
< Lz >= 251. Figure 9b is similar to Figure 9a but for the chaotic orbit of Figure 
8d. Here one can see an asymmetric curve and < Lz >= —27.1. Both orbits were 
calculated for a time period of 100 time units and with n = 10^. 

5. DISCUSSION 

It is well known that the evolution of galaxies is closely associated to the shape 
and the structure of the host halos. Today astronomers believe that dark halos may 
have a variety of shapes and they usually show asymmetries (see, e.g., Dinshaw et 
al. 1998; Oppenheimer et al. 2001; McLin et al. 2002; Penton et al. 2002; Steidel 
et al. 2002). 

In the present paper, the behavior of orbits in a galaxy with an asymmetric dark 
halo component was investigated. In order to keep things simple, we have chosen 
a nearly spherical dark halo component with a small deviation from spherical 
symmetry, introduced by the term —Xx^. The results of this work strongly suggest, 
that small asymmetries in the halo may play an important role on the regular or 
chaotic nature of orbits. 

We have first studied the 2D system using the classical method of the Poincare 
phase plane. Our numerical experiments show that there are three distinct areas 
of chaotic motion when A — 0.015 and = 500: a chaotic layer in the central 
region and two additional chaotic regions, which appear near the unstable periodic 
orbits produced by the secondary resonances. As expected, the LCE was found 



Chaos in a disk galaxy model induced by dark halo 



215 



to be different in each chaotic region. On the other hand, when A = 0.03 and /12 
= 500, the LCE was found to converge to a unique value when the three chaotic 
regions merge to form a large chaotic sea. 

Using the results of the 2D system, we have investigated the 3D potential. We 
have used initial conditions {xcp^o, Zq), where {xo,Pxo) was a point on the phase 
plane. It was found, that all orbits starting in the chaotic regions on the phase 
plane, produce chaotic orbits for all values of zq. Note that in all cases the value 
of Py was found from the energy integral, while y = = 0. Our numerical results 
indicate that the three chaotic components, found in the 2D system, exist also 
in the 3D potential and have also different values of the LCE when A = 0.015. 
Furthermore, when A = 0.03 and the three chaotic components merge, as it was 
indicated in the phase plane of the 2D system, the LCE of the 3D system seems 
to converge to a unique value. This result is in agreement with the outcomes for 
the 3D systems obtained by Cincotta et al. (2006). 

It was also found that orbits starting in the regular parts of the 2D system dis- 
play chaotic motion only when \ zo\ > 4. Orbits with low values oi < >, moving 
near the galactic plane, when approaching the dense and massive nucleus, are scat- 
tered to the halo displaying chaotic motion. This behavior was also observed in 
axially symmetrical disk galaxies with dense massive nuclei (see also Caranicolas 
& Papadopoulos 2003 and references therein). Thus, one can conclude that low 
angular momentum stars are on chaotic orbits in axially symmetrical, as well as 
in triaxial galaxies, hosting massive and dense nuclei. 
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